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Abstract 

We provide a new algorithm for the treatment of deconvolution 
on the sphere which combines the traditional SVD inversion with 
an appropriate thresholding technique in a well chosen new basis. 
We establish upper bounds for the behaviour of our procedure for 
any Lp loss. It is important to emphasize the adaptation properties 
of our procedures with respect to the regularity (sparsity) of the 
object to recover as well as to inhomogeneous smoothness. We also 
perform a numerical study which proves that the procedure shows 
very promising properties in practice as well. 
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1 Introduction 

The spherical deconvolution problem was first proposed by Rooij and Ruym- 
gaart (1991) [H] and subsequently solved in Healy et al. (1998) [3]. Kim 
and Koo (2002) [7| established minimaxity for the L2-rate of convergence. 
The optimal procedures obtained there are using orthogonal series meth- 
ods associated with spherical harmonics. One important problem arising 
with these procedures is their poor local performances due to the fact that 
spherical harmonics are spread all over the sphere. This explains for in- 
stance the fact that although they are optimal in the L2 sense, they cease 
to be optimal for other losses, such as hp losses for instance. 

In our approach we focus on two important points. We aim at a pro- 
cedure of estimation which is efficient from a L2 point of view, as well as 



it performs satisfactorily from a local point of view (for other hp losses for 
instance). 

Deconvolution is an inverse problem and in such a problem there is a 
notable conflict between the inversion part which in presence of noise cre- 
ates an instability reasonably handled by a Singular Value Decomposition 
(SVD) approach and the fact that the SVD basis very rarely is localized 
and capable of representing local features of images, which are especially 
important to recover. Our strategy is to follow the approach started in 
Kerkyacharian et al. (2007) |6] for the Wicksell case , Kerkyacharian et 
al. (2009) [5] for the Radon transform, which utilizes the construction bor- 
rowed from Narkowich Petrushev and Wald (2006) |9j, |8j of a tight frame 
(i.e. a redundant family) staying sufficiently close to the SVD decomposi- 
tion but which enjoys at the same time enough localisation properties to 
be successfully used for statistical estimation (see for instance Baldi et al. 
(2009) Pietrobon et al. (2008) [HISlIOj for other types of applications). The 
construction ^ produces a family of functions which very much resemble 
to wavelets, the needlets. 

To achieve the goals presented above, and especially adaptation to dif- 
ferent regularities and local inhomogeneous smoothness, we essentially use 
a projection method on the needlets (which enables a stable inversion of the 
devonvolution, due to the closeness to the SVD basis) with a fine tuning 
subsequent thresholding process. 

This provides a reasonably simple algorithm with very good perfor- 
mances, both from a theoretical point of view and a numerical point of 
view. In effect, this new algorithm provides a much better spatial adap- 
tation, as well as adaptation to wider classes of regularity. We give here 
upper bounds obtained by the procedure over a large class of Besov spaces 
and any hp losses. 

It is important to notice that especially because we consider different 
hp losses, we provide rates of convergence of new types attained by our 
procedure, which, of course, coincide with the usual ones for L2 losses. 

Again, the problem of choosing appropriated spaces of regularity on the 
sphere in a serious question, and we decided to consider the spaces which 
may be the closest to our natural intuition: those which generalize to the 
sphere case the classical approximation properties of usual regularity spaces 
such as Holder spaces and include at the same time the Sobolev regularity 
spaces used in Kim and Koo (2002) [7]. 

Sphere deconvolution has a vast domain of application; our results are 
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especially motivated by many recent developments in the area of observa- 
tional astrophysics. 

It is a common problem in astrophysics to analyse data sets consisting of 
a number of objects (such as galaxies of a particular type) or of events (such 
as cosmic rays or gamma ray bursts) distributed on the celestial sphere. In 
many cases, such objects trace an underlying probability distribution / on 
the sphere, which itself depends on the physics which governs the production 
of the objects and events. 

The case for instance of ultra high energy cosmic rays (UHECR) iillus- 
trates well the type of applications of our results. Ultra high energy cosmic 
rays are particles of unknown nature which arrive at the earth from appar- 
ently random directions of the sky. They could originate from long-lived 
relic particles from the Big Bang, about 13 billion years old. Alternatively, 
they could be generated by the acceleration of standard particles, such 
as protons, in extremely violent astrophysical phenomena, such as cluster 
shocks. They could also originate from Active Galactic Nuclei (AGN), or 
from neutron stars surrounded by extremely high magnetic fields. 

Hence, in some hypotheses, the underlying probability distribution for 
the directions of incidences of observed UHECRs would be a finite sum of 
point-like sources - or near point like, taken into account the deflection of 
the cosmic rays by magnetic fields. In other hypotheses, the distribution 
could be uniform, or smooth and correlated with the local distribution of 
matter in the universe. The distribution could also be a superposition of the 
above. Identifying between these hypotheses is of primordial importance 
for understanding the origin and mechanism of production of UHECRs. 

Of course, the observations of these events (Xj's in the sequel) are 
always most often perturbated by a secondary noise (e,) which leads to 
the deconvolution problem described below. Following Healy et al. (1998) 
[3], Kim and Koo (2002) [7j, the spherical deconvolution problem can be 
described as follows. Consider the situation where we observe Zi, . . . , 
N i.i.d. observation with 

Zi = EiXi, (1) 

where the e^'s are i.i.d. random elements in 5*0(3) (the group of 3 x3 
rotation matrices) , and the Zj's and Xj's are i.i.d. random elements of 
§^ (two-dimensional unit sphere of M^) random elements, with Si and Xi 
assumed to be independent. We suppose that the distribution of resp. X, 
Z and e are absolutely continuous with the Haar measure of resp. §^ 
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and 50(3) with the densities fz, fe, fx- 
Then, 

fz = fe*fx, (2) 

where * denotes convolution and is defined below. In the sequel, fx will be 
denoted by / to emphasize the fact that it is the object to recover. 

The following paragraph recall the necessary definitions. It is largely 
inspired by Kim and Koo (2002) [7] and and Healy et al. (1998) [3]. 

2 Some preliminaries about harmonic anal- 
ysis on SO{3) and 

We will provide a brief overview of Fourier analysis on 5*0(3) and Most 
of the material can be found in expanded form in Vilenkin (1969) [15], 
Talman (1968) [12], Terras (1985) [13], Kim and Koo (2002) [7], and Healy 
et al. (1998) [3]. Let 

(cos — sin \ / cos 9 sin 9 \ 

sin0 COS0 , a{9) = 10, 
1/ \ -sin^ cos9 J 

where , G [0,27r), 9 G [0,7r). It is well known that any rotation matrix 
can be decomposed as a product of three elemental rotations, one around 
the z-axis first with angle ip, followed by a rotation around the y-axis with 
angle 9, and finally another rotation again around the z-axis with angle 0. 
Indeed, the well known Euler angle decomposition says that any g G 50(3) 
can almost surely be uniquely represented by three angles (0, 9, '0),with the 
following formula (see Healy et al. (1998) [3] for details) : 

g = u{<f))a{9)u{^), (3) 

where G [0,27r), 9 G [0,7r), tp G [0,27r). Consider the functions, known 
as the rotational harmonics, 

DU<P,0,,P) = e-(-'^+"'^)pl(cos^), (4) 

where the generalized Legendre associated functions for —l<m, n < 
/, / = 0, 1,... are fully described in Vilenkin (1969) [15]. The functions 
D^^^ ioT —I < m, n < I, 1 = 0, 1, . . . are the eigenfunctions of the Laplace 
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Beltrami operator on 50(3), hence, a/2/ + l-D^n; —l^m, n < I, I = 
0, 1, ... is a complete orthonormal basis for L2(S'0(3)) with respect to the 
probability Haar measure. In addition, if we define the {21 + 1) x {21 + 1) 
matrices by 

D\g) = [D^M], (5) 

where for — / < m, n < I, 1 = 0, 1, . . . and g G S0{3), they constitute the 
collection of inequivalent irreducible representations of 50(3) (for further 
details see Vilenkin (1969) [I5]). 

Hence, for / G L2(S'0(3)), we define the rotational Fourier transform 
on 50(3) by 



fLn= / fi9)Dlni9)d9, (6) 

JS0{3) 

where again we think of ^ as the matrix entries of the {21 + 1) x {21 + 1) 
matrix 

/ = [fmn]~l<m, n<U / = 0, 1, . . . 

and dg is the probability Haar measure on 50(3). The rotational inversion 
can be obtained by 



I —l<m, n<l 

= E E tuDLig-') (7) 



I —Km, n<l 



([7]) is to be understood in L2-sense although with additional smoothness 
conditions, it can hold pointwise. 

A parallel spherical Fourier analysis is available on Any point on 
§^ can be represented by 

LJ = {cos(f) sin 9, sin sin 9, cos 6*)*, 

with , G [0, 2n), 9 G [0, n). We also define the functions : 



>1H = yLid,<t>) = {-lr^r^^^^^^PL{cos9)e 



Att {I + m) 



for -/ < m < /, / = 0, 1, . . ., G [0,27r), 9 G [0, vr) and where P4(cos^) 
are the associated Legendre functions. The functions obey 

yUo,^) = {-irYi{9,ct>). (9) 
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The set {Yj^, —l<m < I, I = 0, 1, . . .} is forming an orthonormal basis 
of L2(§^), generally referred to as the spherical harmonic basis. 

Again, as above, for / G L2(§^), we define the spherical Fourier trans- 
form on by 

t = [ fiu;)nMdu^, (10) 

where du is the probability Haar measure on the sphere S^. The spher- 
ical inversion can be obtained by 

/H = E E f'mYii^)- (11) 

I -l<rn<l 

The bases detailed above are important because they realize a singular value 
decomposition of the convolution operator created by our model. In effect, 
we define for G L2(5'0(3)), / G L2(§^) the convolution by the following 
formula: 

fe*fi^)= fe{u)f{u-^Uj)du 
J 30(3) 

and we have for all — / < m < /, / = 0, 1, . . ., 

ifehm- (12) 



(/e * f )m ^ ^ fe,mnfn 



2.1 The SVD Method 

The singular value method (see Healy et al. (1998) [3J and Kim Koo (2002) 
[7]) consists in expanding / in the spherical harmonics basis Yj^ and esti- 



mating the spherical Fourier coefficients using the formula above (12). We 
get the following estimator of the spherical Fourier transform of /: 

^ N I 

/m^ •= E E f£'^,mrXni^j) (13) 



N 



provided, of course, that these inverse matrices exist, and then the estimator 
of the distribution / is 



N 



/"'(^) = E E f'rfyii-^), (14) 



1=0 m=~l 
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where depending on the number of observations has to be properly se- 
lected. 

3 Needlet construction 

This construction is due to Narcowich et al. (2006) |B]. Its aim is essen- 
tially to build a very well localized tight frame constructed using spherical 
harmonics, as discussed below. It was recently extended to more general 
euclidean settings with fruitful statistical applications (see Kerkyacharian 
et al. (2007) [6]), Baldi et al. (2009) Pietrobon et al. (2008) [IlilllO]. As 
described above, we have the following decomposition: 

oo 

L2(c/x) = 0Hz, (15) 

1=0 

where is the space of spherical harmonics of S^, of degree / (which di- 
mension is 2Z + 1). 

The orthogonal projector on EI; can be written using the following the 
kernel operator 

V/ eh^idx), PuJ{x)= [ Li{{x,y))f{y)dy (16) 

where {x,y) is the standard scalar product of M.^, and Li is the Gegenbauer 
polynomial with parameter | of degree /, defined on [—1,-1-1] and normalized 
so that 

j ^L,{t)L,{t) = ^ (17) 

Let us point out the following reproducing property of the projection 
operators: 

/ Li{{x,y))Lk{{y,z))dy = 6i^kLi{{x,z)) . (18) 

The following construction is based on two fundamental steps: Littlewood- 
Paley decomposition and discretization, which are summarized in the two 
following subsections. 
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3.1 Littlewood-Paley decomposition 

Let be a function on M, symmetric and decreasing on M"*" supported 
in 1^1 < 1, such that 1 > 0(0 > and 0(0 = 1 if |C| < |- 

&'(O = 0(|)-0(O>o 

so that 

i>o 

Remark that 6(0 7^ only if | < |0 < 2. Let us now define the 
operator 

Aj — ^;>o and the associated kernel 

We obviously have: 

J 

V/ e L2(S2), / = hm Lo(/) + J] A,(/) ■ (20) 

■/— »oo ' 

and if Mj{x,y) = T.i>oKij)Li{{x,y)), then 

Aj{x, y) = j Mj{x, z)Mj{z, y) dz . (21) 

3.2 Discretization and localization properties 

Let us define 

I 

m=0 

the space of the restrictions to of the polynomials of degree less than 

I. 

The following quadrature formula is true: for alH e N there exists a 
finite subset of S"^ and positive real numbers > 0, indexed by the 
elements of such that 
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V/G^z, / fix)dx=J2\fiv)- (22) 
Then the operator Mj defined in the subsection above is such that: 

Z ^ Mj{x,z) e ^[2.+i] , 

so that 

z ^ Mj{x,z)Mj{z,y) G ^[2^+2] , 
and we can write: 

A,{x,y) = j M,{x,z)M,{z,y)dz= ^ X^Mj{x,r])Mj{r],y) . 
This imphes: 

Ajf{x) = j Aj{x,y)f{y)dy = j Yl \M,{x,r])M,{'n,y)f{y)dy 



= \^vMjix,r]) J ^^M,{y,ri)f{y)dy . 

We denote 

<^[2J+2] = ^j, '^j,vi^) ■= Mj{x,ri) for rj e 3fj , 

It can also be proved that the set of cubature points can be chosen 
so that: 

l2^^ <i^3fj <c2'^^ (23) 
for some c > 0. Actually in the simulations of ^we make use of some 



sets of cubature points such that = 12.2^-' exactly. It holds, using (20) 
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The main result of Narcowich et al. (2006) is the following local- 
ization property of the ipj^rj, called needlets: for any k E N there exists a 
constant Ck such that, for every ^ G 

where d is the natural geodesic distance on the sphere {d{^, rj) = arccos(?7, ^) 
). In other words needlets are almost exponentially localized around their 
associated cubature point, which motivates their name. 

A major consequence of this localization property can be summarized 
in the following properties which will play an essential role in the sequel. 

For any 1 < p < oo, there exist positive constants Cp, Cp, c, C and Dp 
such that 



||$:A,^,,|[<c2^^-(i-^)($:|A, 

II Yl 

II Urjlpjri 



< 



DpJ2\>^ 



I/tt 



< C sup \un\2^ 



(25) 
(26) 

(27) 
(28) 
(29) 



To conclude this section, let us give a graphic representation of a spher- 
ical needlet in the spherical coordinates in order to illustrate the above 
theory. In the following graphic, we chose j = 3 and rj = 250. 
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Figure 1: A spherical needlet. 
3.3 Besov spaces on the sphere 

The problem of choosing appropriated spaces of regularity on the sphere in 
a serious question, and we decided to consider the spaces which may be the 
closest to our natural intuition: those which generalize to the sphere case 
the classical approximation properties used to define for instance Sobolev 
spaces. In this section we summarize the main properties of Besov spaces 
which will be used in the sequel, as established in |8]. 
Let / : §^ ^ M be a measurable function. We define 

E,{f,n)= mi II/-P1U, 

the distance between / and the space of polynomials of degree k. 
The Besov space B^ ,^ is defined as the space of functions such that 

feu and (j2ik'Ek{f,7r)yj:) ' < +00. 

fc=0 

Remarking that k — * Ek{f, vr) is decreasing, by a standard condensation 
argument this is equivalent to 
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/ e and (^(2^'^E2.(/,7r))'-) ' < +00. 
i=o 

and the following theorem states that as it is the case for Besov spaces 
in R'^, the needlet coefficients are good indicators of the regularity and in 
fact Besov spaces of are Besov bodies, when expressed using the needlet 
expansion. 

Theorem 1. Let 1 < tt < +00, s > 0, < r < +00. Let f a measurable 
function and define 

provided the integrals exist. Then f belongs to B^ ,^ if and only if for 
every j = 1, 2, . . . , 

where {6j)j e ir- 

As has been seen above, 

c22^'(^-^) < < C22^'(5-i), 

for some positive constants c, C, the Besov space 5^ turns out to be 
a Banach space associate to the norm 

WfWns^^ := ||(2^[^+'(^-^)l||(/3,,),e^,|kJ,>o|k. < 00, and (30) 

Using standard arguments (reducing to comparisons of Iq norms) , it is easy 
to prove the following embeddings: 

C B;^, forp<7r 

B' C Bp'r^^^'^^ for TT < p and s > 2(- - -). (31) 

Moreover, it is also true that for s > f , if / belongs to 5*^, then it is 
continuous, and as a consequence bounded. 

In the sequel we shall denote by B^ ,^{M) the ball of radius M of the 
Besov space B^^.. 
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4 Needlet algorithm: thresholding needlet 
coefficients 



The first step is to construct a needlet system (frame) {ipjn : r] & > 
— 1} as described in section 3. 

Tlie needlet decomposition of any / G L2(§^) takes the form 

j r)(^S] 

Using Parseval's identity, we have (3jy^ = (/, i^jr^)h2(s^) = J2im fLi^%' with 
= and = (^,„Fj,). 

Thus 

Piv = E /^""^J^' (32) 



Im 



is an unbiased estimate of We recall that has been defined in (13). 

Notice that from the needlet construction (see the previous section) 
it follows that the sum above is finite. More precisely, ip'j^ 7^ only for 
2^-1 < / < 2^+\ 

Let us consider the following estimate of /: 

J 

where t is a thresholding operator defined by 

t{l3jr^) = Pjr,H\Pjv\ ^ f^^NWjl} with (33) 

-l = ^'EiE^5^/'-wP- (35) 

In m 

Here k is a tuning parameter of the method which will be properly 
selected later on. M is such that ||/||oo < M. Notice that the thresholding 
depends on the resolution level j through the constant Uj which will also 
be specified later on, and the same with regard to the upper level of details 
J. 
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4.1 Performances of the procedure 

Theorem 2. Let 1 < p < oo, u > 0, and let us assume that 



In m 



Let us take > 16p and 2"^ = d[tN\ *''+^' with t^v m (34) anc? d is a 
positive constant. 
Then 

^/tt > 1, s > 2/7r, r > 1 (wi/i the restriction r < ix if s = (z^+l)(- — 1)), 
t/iere exists a constant C such that: 

sup E||/ - fWl < Ci\ogiN)r-'[N-'/'Vhi(N)r, (37) 



where 



= — z/.>(z. + l)(^-l) 

S + + 1 TT 

s-2/7r + 2/p .2 , ^, 

= ^ o/ ii > if - <s<{u + l){^-l). 

S + U — 2/tT + 1 IT 71 

The proof of this theorem is given in section |6l 



Remarks 



1. The rates of convergence found here are standard in inverse problems. 
They can be related to rates found in Kim and Koo (2002) in the same 
deconvolution problem, with a L2 loss and constraints on the spaces 
comparable to Bg2{M). In the deconvolution problem on the interval, 
similar rates are found even for hp losses (with standard modifications 
since the dimension here is 2 instead of 1): see for instance Johnstone 
et al. (2004) [4]. These results are proved to be minimax (see Kim 
and Koo (2002) |7j) up to logarithmic factors, for the case p = 2 with 
a Bg2{M) constraint on the object to estimate. 

2. It is worthwhile to notice that the procedure is adaptive, meaning that 
it does not require a priori knowledge on the regularity (or sparsity) 
of the function. It also adapt to non homogeneous smoothness of 
the function. The logarithmic factor is a standard price to pay for 
adaptation. 
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3. The parameter v appearing here is often called degree of ill-posedness 



of the problem (DIP). It appears here through condition (36) which 
is essential in this problem. In [7] for instance, and very often in 
diverse inverse problems, this DIP parameter is introduced with the 
help of the eigenvalues of the operator (i.e. here the discrepancy of 
the coefficients of in its expansion along the spherical harmonics). 



In the following subsection, we prove that (36) is in fact a consequence 



of the standard 'ordinary smooth' condition. 



4.2 Condition (36) and the smoothness of /g 



Following Kim and Koo (2002) [7] (condition 3.6), we can define the smooth- 
ness of spectrally. We place ourselves in the 'ordinary smooth' case 

\\{'!\Y^\\ov<d^V and \\'!%p<dr'' as / ^ cx), (38) 

for some positive constants (io, d\ and nonnegative constant and where 
the operator norm of the rotational Fourier transform j\ is defined as 

llfll = sup 

WleWov ||T|| 5 

£,[ being the (2Z + l)-dimensional vector space spanned by : —I < m < 
I}. 



The following proposition states that condition (35) is satisfied in the ordi- 
nary smooth case by the needlets system. 

Proposition 1. // H/fHop < dol''^ , then 

In m 

Proof. Since i)]";; ^ only for 2^-^ < I < 2^+1, 

2j+l I 2J+1 



1=23-^ n=— I m 1=21' 



2J+i 
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which proves the result using inequahty (25). 



We now give a brief review of some examples of smooth distributions which 
are discussed in depth in Healy et al (1998) |3] and Kim and Koo (2002) 

4.2.1 Rotational Laplace distribution 

This distribution can be viewed as an exact analogy on SO (3) of the Laplace 
distribution on M. Spectrally, for some > 0, this distribution is charac- 
terized by 

fUn = a+p'Ki + i)r'5^n, (39) 

for —l<m,n<l and / = 0, 1, and where 6mn = 1 if m = n and 
otherwise. 



4.2.2 The Rosenthal distribution 

This distribution has its origin in random walks in groups (for details see 
Rosenthal (1994) [llj). 

If one considers the situation where is a p-fold convolution product 
of conjugate invariant random for a fixed axis, then, Rosenthal (1994) [TT] 
(p. 407) showed that 

^ ( sin(/ + l/2)g y 
for —l<m,n<l and Z = 0, 1, ... and where < 6 < 7i and p > 0. 



5 Practical performances 

In this section we produce the results of numerical experiments on the 
sphere The sets of cubature points ^jrj in the simulations that follow 
have been generated with the HEALPix pixelisation. For each resolution 
level j, the HEALPix pixelisation gives 12. 2^-^ cubature points. 
In the two examples below we considered samples of cardinality = 1500. 
The maximal resolution level J is taken such that J = (l/2)log2 (i~/v)' 
In order not to have more cubature points than observations we set J = 3 
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for = 1500. We recall the expression of the estimate of the needlets 
coefficients of the density of interest: 

I I N 

i=2J-i m=-l n=-l u=l 

where the cubature weight is equal to A^^ = 47r/(12.2^-'). 
We replace the rotational Fourier transform (/^)mn := flmn (defined in (6)) 
by its empirical version. Indeed, we suppose that the noise is unknown, 
this is a situation which is very likely to occur for instance in the context 
of astrophysics. 

We precise again that denotes the (m, n) element of the matrix 

{fl)~^ := /^-i which is the inverse of the (2/ + 1) x (2/ + 1) matrix (fl). 
In order to get the empirical version of /^-i we have first to 

compute the empirical matrix (Z^'^) then to inverse it to get the matrix 
ife^)~^ ■~ /e^- '^^^ {m,n) entry of the matrix (/J'^) is given by the 
formula: 

N 

fe,mn N Z)^^„(£j), 

i=i 

where the rotational harmonics -D^„ have been defined in Q. The e^-'s are 
i.i.d realizations of the variable e G 5*0(3). 

For the generation of the random variable e G SO (3), we chose the Oz as 
the rotation axis and an angle following a uniform law. Various cases 
are investigated for the angle 0, that is to say a uniform law with different 
supports such as [0,7r/8], [0,7r/4], [0,7r/2]. The more the support of the 
uniform law is large, the more the level noise will be intense. 
This particular choice of rotation matrix entails in the decomposition of an 
element of SO (3) (see formula ([3])) the angles 'ip and 6 to be both equal to 
zero. For this specific setting of perturbation, we deduce the following form 
for the rotational harmonics: 

where ~ t/[0, a] and a is a positive constant which will be specified later. 
For the reconstruction of the density /, we recall that 



^ J 12.22J 

f |S2| + E E '^i'?^i'?-^{|/3j>,|>«:tn(Tj}' 



j=0 ri=l 
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the expression of aj is given in (35), and t]y = ^y\og{N) / N . The constant 
M which appears in the expression of aj is such that ||/||oo < M. 



Example 1. In this first example, we consider the case of the uniform 
density f = It is easy to verify that Pjr, = {f,4'jrf)L'2 = for every j 
and every t]. Following Baldi et al. (2009) [Ij, a simple way of assessing 
the performance of the procedure is to count the number of coefficients 
surviving thresholding. 
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Table 1: number of coefficients surviving thresholding for various values of 
Ko, 0~f/[O,7r/8]. 
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Table 2: number of coefficients surviving thresholding for various values of 

Kq, (f) ~ U[0,Tt]. 

We precise that both in the cases of an angle following a law f/[0, 71/8] 
or f/[0,7r], all the coefficients are killed for kq = 0.38. Accordingly, we can 
conclude that the thresholding procedure based on spherical needlets is very 
efficient. 

Example 2. We will now deal with the example of a density of the 
form f{Lu) = ce-^l'^-^il", with = (0,1,0) and c = 1/0.7854. With this 
choice of density, it turns out that ||/||oo = 1/0.7854 = 1.2732. Hence we 
can set M such that M = 1.2732. The graph of / in the spherical coor- 
dinates ($,e) ($ = longitude, < $ < 27r, 6 =colatitude, < 9 < vr) 
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is given in Fig 2. We also plot the noisy observations for different cases 
of perturbations. We remark that for rather big rotation angles such that 
(f) ~ C/[0,7r/4] or ~ C/[0, 7r/2], the observations tend to be spread over a 
large region on the sphere and not being concentrated in a specific region 
any more. Consequently, denoising might prove to be difficult. In the con- 
text of the deconvolution on the sphere, a large amount of noise corresponds 
to a rotation around the Oz axis with a large angle. 

A fundamental issue in the field of astrophysics is to detect the place 
of the peak of the bell which in the present density case is localized in 
(0 = 7r/2, $ = 7r/2). For each case of noise, we plot the observations both 
on the sphere and on the fiattened sphere and give the reconstructed den- 
sity in the spherical coordinates. 
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Figure 3: Observations (p ~ U[0,tt/8]. 




Figure 4: Observations (p ~ f/[0,7r/8]. 
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Figure 5: The estimated density, band with kq = 0.43. 




-1 



Figure 6: Observations ~ f/[0,7r/4]. 
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Figure 7: Observations ~ [/[0,7r/4]. 
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Figure 11: The estimated density, bandwith kq = 0.56. 

At a closer inspection, we notice that the position of the peak of the 
estimated bell is pretty well localized whatever the amount of noise. Never- 
theless in the case of the law U[0, vr/2], the longitude coordinate of the peak 
tends to slightly move, its colatitude coordinate remaining well localized. 
Therefore even if in the case of rather big rotations such that (p ~ f/[0, 7r/4] 
and (p ~ f/[0,7r/2], our estimation procedure allows us to detect the posi- 
tion of the peak. Of course, one remarks that the base of the bell tends to 
become a bit larger when the noise increases, this is due to the fact that 
the observations are not concentrated in a specific region any longer, but 
the genuine form of the density is well preserved. 



6 Proof of Theorem [2 

In this proof, C will denote an absolute constant which may change from 
one line to the other. 

We begin with the following proposition 

Proposition 2. For any p > 1, 

m,-Pjv\' <c[|]' (41) 

n\Pj^ - Pjrjl > cr.Kt^] < ciV-^W (42) 

= -^2/4 (43) 
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To prove (41), we will use Rosenthal inequality: 

N 



f3j„ (3. 



X) 



Hence, 

VarG{Z) < I \GiLu)\'' fz{u;)dLU < M I \G{uj)\''duj, 



since ||/z||oo < M (this is a consequence of the fact that s > Ijix which in 
turns implies that / is continuous, and so is fz on the sphere). Then using 
Parseval: 



I \G{yj)\^d^ = Y^\Y^^fju 

-^^^ In m 



|2 _ 2 



Moreover, using Cauchy Schwarz inequality, and noticing again that ^ 
only for 2^-^ < I < 2^+^: 

In m Ini 

since J2m l^mi^)]'^ = because of the addition formula ( see Kim and 
Koo (2002) |7]) inequality (6.1)). Then using Rosenthal inequality, we get 
(we recall that 2^ < iV^/^ j < J): 

1 ^ 1 1 

nj^Y.^{Z,)-^G{Z,)Y < C[—a^{a,2^/'r-'+ — {Na^)^ < Ca]N~^l\ 

1=1 



To prove (42), we use Bernstein inequality: 

This ends up the proof of the proposition. 

Now, to get the result of Theorem [2| we begin by the following decom- 
position : 

J 

ni -m < 2^"Heii e e m,) - /3,.)^,.ii^ + ii E E f^Mi;} 

=: / + // 
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The term // is easy to analyse, as follows. We observe first that since 
S^^,(M) C 5^,,(M ) for vr > p, this case will be assimilated to the case 
TT = p and from now on, we will only consider tt < p. Since / belongs to 



B^j.{M), using the embedding results recalled above in (31 ), we have that / 



also belongs to Bp^r " {M'), for some constant M' and for tt < p. Hence 



j>J rieS'j 



-J[s-2(i-i)] 



s-2{---) 

Then we only need to verify that — ^ is always larger that /x, which is 
not difficult. 

Indeed, on the first zone s > l)(p/7r — 1). So, s + u+l > (z/ + l)2 which 
entails that , < , J,.p ■ We need to check that s — 2(- — -) > —. We 

have that s - - + - - - = 2(f - - > since s > - and » > vr. 

■K p p ^2 ' ^-K pi TT 

On the second zone, we obviously have that — ^'^^ ^ is always larger that 

^ s-2/7r+2/p 
s+i/-2/7r+l' 

Bounding the term I is more involved. Using the triangular inequality 



together with Holder inequality, and property (28) for the second line, we 
get 

J 

i=-l r^eiTj 
J 

i=-ir,eir, 
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Now, we separate four cases : 



+I{\(3jr^\ < KtNCTj} 



HlPrril > pN(rj} + I{\Pj^\ < pN(rj} 

+ \Pjv\''\\^Jv\\lI{\P]v\ < K^tN(Tj}\l{\Pj^\ > 2KtNa,} 



< 2KtNaj} 

= : Bb + Bs + Sb + Ss. 

Using Proposition [2} we have 



And, using in addition (25) and (36) 



J 



j=-lrjG3fj 

J 



Now, since / belongs to B^ j.{M), using again Proposition |2] and (31), we 
have 
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i=-lr;eJ'j- 



12N-^ I{\P,r,\>2KtM(rj} 



i=-i 

It is easy to check that in any cases ii k > 1 — p the terms Bs and Sb 
are smaller than the rates announced in the theorem. 

We have using Proposition |2| (25) and condition (36) for any p > -2 > 0: 



i=-i 
J 



Also, for any j9 > 2; > 



\p-z 



i=-i V&3fj 

So in both cases we have the same bound to investigate. We will write 
this bound on the following form (forgetting the constant) : 



JO 



A+B -.^ t„'-"[ 



i=jo+i 
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The constants Zi and jo will be chosen depending on the cases, with the 
only constraint p > -Zj > 0. 

We recall that we only need to investigate the case p > tt, since when 
p < vr, Bl,{M) C B;,{M'). 

Let us first consider the case where s>(i/ + l)(^ — 1), put 



S + + 1' 



and observe that on the considered domain, q <7T and p > q. In the sequel 
it will be useful to observe that we have s — {u + 1)(| — 1). Now, taking 
Z2 — TT, we get: 

J 



Now, as 



and 



h vi = s + l 

q 7^ q TT 



with (rj)j e (this last thing is a consequence of the fact that / e 
B^ ,^{M)), we can write : 

i=io+i 

< Ct^P-'^2^°^(^"?)(^+^\ 

The last inequality is true for any r > 1 if n > q and for r < tt if tt = g. 
Notice that vr = g is equivalent to s = (z/ + 1)(^ — 1). Now if we choose jo 
such that 2-'°9'^^+^) t^'^ we get the bound 

which exactly gives the rate announced in the theorem for this case. As 
for the first part of the sum (before jo); we have, taking now Zi = q, 
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with q < 71 (and also q < p since we investigate the case p > tt), so that 
[w E^eiT,. < iw E^e^S- \Pjv\"]^^ we get 



JO 

-1 ryeiTj 

io 

-1 »?e^3 

io 



< 



The last two lines are valid if q is chosen strictly smaller than q (this is 
possible since n > q). 

Let us now consider the case where s<(z/ + l)(2 — 1), and choose now 

Z/ + 1 - ^ 

p 

q = p TT-^ — . 

In such a way that we easily verify that p — q = p ^i^J]^^^lj^ , q — tt — 
ip-^)ii+u)-ns Q Furthermore we also have s + 1 - ^ = 2 - ? + ^^(2 _ i). 

Hence taking zx — 'k and using again the fact that / belongs to ^.(M), 

io 

-1 ri^Sfj 

io 
-1 



IS IS 11 ue siiiue -r J- — ^ 

2 



This is true since v + 1 — - is also strictly positive since v + 1 > > 
- > -. If we now take 2^°9^^'^^~i^ ~ tN'^ we get the bound 



P 
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which is the rate announced in the theorem for this case. 
Again, for we have, taking now Z2 = q > q{> n) 

J 

j=jo+i n&sfj 

But 

and s + 1 - - ^ ^ - - + - 1), hence 

i=io+i 
< W, 

which completes the proof of Theorem 2. 
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